Patterns and drivers of soil surface‐dwelling Oribatida diversity along an altitudinal gradient on the Changbai Mountain, China

Abstract Distribution patterns of biodiversity and environmental interactions are dominant themes in ecology. In montane ecosystems, biodiversity is closely associated with altitudinal gradients. However, studies of biodiversity in montane ecosystems are focused on plants and vertebrates, with relatively less on invertebrates. Here, the present study used a Vortis arthropod suction sampler to explore the biodiversity patterns of soil surface‐dwelling Oribatida and their drivers along an altitudinal gradient (600, 800, 1600, 2000, and 2300 m) from typical temperate forests, evergreen coniferous forests, subalpine birch forests to alpine tundra on the north slope of Changbai Mountain, Northeast China. Trichoribates berlesei, Platynothrus peltifer, and Oribatula tibialis were the dominant soil surface‐dwelling species on Changbai Mountain. Generally, alpha diversity and beta diversity of soil surface‐dwelling Oribatida decreased with the rising altitude, with a peaking density value at 2000 m. The result of beta diversity showed that the structures of community were more influenced by the species turnover component than the nestedness component. Nonmetric multidimensional scaling (NMDS) ordination showed that the community structure of soil surface‐dwelling Oribatida varied significantly along the altitudinal gradient. The variance partitioning showed that the elevation and climatic conditions determined the soil surface‐dwelling Oribatida community. Spatial filtering represented by geographic and elevation distances was particularly associated with soil surface‐dwelling Oribatida community variation between altitudes on Changbai Mountain. However, the variation of the Oribatida community between adjacent altitudes was only associated with geographic distance. Our study provides supportive evidence for the biodiversity analyzing of soil surface‐dwelling Oribatida in montane ecosystems along an altitudinal gradient.


| INTRODUC TI ON
Montane ecosystems comprise approximately 25% of the biodiversity of all terrestrial ecosystems (Miller & Spoolman, 2011;Rahbek et al., 2019). The most significant difference between montane and other ecosystems is the divergence in vertical zonation along the altitudinal gradient. The altitudinal gradient can be regarded similar to the climatic gradient, compressing the numerical variety of important factors such as temperature, precipitation, ultraviolet radiation, air pressure, wind velocity, and habitat characteristics into a relatively short geographical distance (Barry, 1992;Mayor et al., 2017;Reusch & Wood, 2007;Sundqvist et al., 2013). The altitudinal gradient significantly affects montane base resource availability and turnover efficiency (Gentili et al., 2013;Odland, 2009). Montane ecosystems provide a unique opportunity to study patterns and drivers of biodiversity along the altitudinal gradient.
Currently, studies of biodiversity patterns along the altitudinal gradient focus mainly on alpha diversity and less on beta diversity (Li et al., 2019;Winkler et al., 2018). Beta diversity can reflect the variation of biological communities along the habitat gradient (Bishop et al., 2015). Species turnover and nestedness are fundamental components of spatial community dissimilarity to measure the difference in species composition, species turnover replaces an existing species with a different species at a new site, and nestedness is the extent of species-poor assemblages that constitute subsets of species-rich assemblages (García-Navas et al., 2020;Koleff et al., 2003;Wang et al., 2018). Thus, combining alpha diversity and beta diversity will support a better understanding of the biological distribution patterns of montane ecosystems, as well as predict the variation of the biological community to a changing environment (Gaston, 2000;Sundqvist et al., 2013).
Because of the workload involved in sampling invertebrates of montane ecosystems is enormous and complex, studies are often less focused on the variation pattern of invertebrates as compared to those for plants and vertebrates (Freeman, 2020;Sharma et al., 2019). Surface litter layer is a critical component of ecosystems and researchers have increasingly recognized the importance of the surface litter layer in recent years. The physical factors at the soil surface above the litter layers show more fluctuations than the soil environment (Otto & Svensson, 1982). Compared to soil invertebrates, soil surface-dwelling invertebrates are closely related to plant diversity, and the variety of litter provides suitable microhabitats for soil surface-dwelling invertebrates (Hasegawa et al., 2006;Xu et al., 2020). The high rate of litter turnover along an altitudinal gradient affects the exposure risk and food availability for invertebrates (Olson, 1994;Valencia-Cuevas et al., 2020).
Oribatida are the most abundant microarthropods in forest ecosystems (Illig et al., 2010). They can directly accelerate organic matter decomposition, nutrient cycling through digestion and decomposing litter and play a critical role in soil microstructure formation (Bradford et al., 2002;González & Seastedt, 2001;Huhta, 2007;Lussenhop & Bassirirad, 2005;Wu et al., 2015;Yin et al., 2019). The diversity of Oribatida is influenced by a long-term combination of factors (Mumladze et al., 2015). The inconsistent patterns of Oribatida along an altitudinal gradient are subject to environmental factors, geographical distance, and specific species characteristics (Illig et al., 2010;Marian et al., 2018). Influenced by these factors, the altitudinal diversity patterns of Oribatida generally decline linearly (Mumladze et al., 2017) or peak at an intermediate altitude, where the primary productivity is the highest .
The Changbai Mountain with its peak of 2691 m is the tallest in Northeast China with a significant vertical zonation. The vertical vegetation distribution on Changbai reflects the horizontal vegetation from temperate to boreal zones in Eurasia (Shao, 2011). This study addresses two questions, including: (1)  To study these research questions, we investigated soil surfacedwelling Oribatida from 600 to 2300 m on Changbai Mountain. We hypothesize that (1) the distribution pattern of soil surface-dwelling Oribatida is affected by deterministic and stochastic factors (Mumladze et al., 2015). Therefore, alpha diversity and beta diversity of soil surface-dwelling Oribatida community would decrease with rising altitude, with a peak of diversity at mid-altitudes. (2) Terrestrial invertebrates typically have a narrow altitudinal range of distribution (Brühl et al., 1999). Therefore, species turnover will be the main process of shaping the community structure of soil surface-dwelling Oribatida along the altitudinal gradient. (3) The process of community construction related to environmental filtering and spatial filtering is essential to the dynamics of species distribution patterns (Heino et al., 2015;Peters et al., 2016), montane ecosystems provide variable environmental gradients, geographic distances, and elevation distances, we hypothesize that the variation in soil surfacedwelling Oribatida community in different altitudes would be more related to spatial filtering and climate change driven by altitude variation rather than other environmental factors.

| Study area
This study was conducted on the Changbai Mountain Nature Northeast China. The region is one of the few well-protected natural ecosystems (Shen et al., 2013). The region is characterized by a temperate continental mountain climate with a mean annual temperature of 3.6°C. The mean annual precipitation is 700 mm, with more than 80% of the precipitation occurring between June and September (Chen et al., 2011;Wu et al., 2002;Yu et al., 2013).

| Sampling methodology
We investigated soil surface-dwelling Oribatida assemblages at five altitudes: 600 m (a mixed coniferous and broad-leaved forest), 800 m The sampler can filter extraneous materials, shorten the sampling periods, and is a standard suction sampler in scientific research for invertebrates such as Collembola, Acari, Araneae, and others. (Bohan et al., 2011;Fitzgerald & Solomon, 2004;Stewart, 2002;Zentane et al., 2016). The soil surface-dwelling fauna was preserved in 95% alcohol and transported to the laboratory for identification. After submerging the Oribatida into lactic acid to make them fade, we identified and counted them under a stereo stereomicroscope with 6.3 to 80× magnification (SMZ1270, Nikon) (Niedbała & Liu, 2018;Ryabinin et al., 2018;Yin, 1998).
The coordinates of latitude, longitude, and elevation were re-  Table 1. We used WorldClim version 2 to download the climate data (Fick & Hijmans, 2017) and extract the mean annual temperature (MAT) and mean annual precipitation (MAP) data by latitude and longitude coordinates.
F I G U R E 1 Geographic locations of the sampling sites along an altitudinal gradient on the Changbai Mountain, Northeast China.

| Statistical analysis
Based on the sampling radius of the Vortis arthropod suction sampler and species abundance, the species density of 1 m 2 was calculated.
Based on the abundance data, the soil surface-dwelling Oribatida were classified into dominant species (>10% of the total number), common species (1-10% of the total number), and rare species (<1% of the total number) (Tan et al., 2010).
To explore altitudinal alpha diversity patterns of Oribatida assemblages, richness, Shannon-Weiner index (H′), and Simpson index (λ) were calculated (Shannon, 1948;Simpson, 1949). After testing the normality of the data, we chose the Kruskal-Wallis test to detect the difference in soil surface-dwelling Oribatida diversity between different altitudes. Multiple linear regression was used to analyze the relationship between the density and alpha diversity of soil surface-dwelling Oribatida community and altitude. Alpha diversity and multiple linear regression were calculated using the "diversity" and "lm" functions. The alpha diversity was calculated with the "vegan" package (Oksanen et al., 2013).
To assess differences in the Oribatida assemblage beta diversity pattern along the altitudinal gradient, we calculated pairwise beta diversity using the Sørensen dissimilarity index. Three matrices (Sørensen pairwise dissimilarity, β sor ; Simpson pairwise dissimilarity, β sim ; Nestedness-resultant dissimilarity, β sne ) of each altitude were constructed. The β sor is overall beta diversity, β sim and β sne are turnover and nestedness components, respectively, and pairwise beta diversity was calculated with the "beta.pair" function in the "betapart" package (Baselga, 2010).
To visualize the differences between soil surface-dwelling Oribatida community compositions along the altitudinal gradient, we calculated the non-metric multidimensional scaling (NMDS) ordination based on Bray-Curtis dissimilarity matrix and analysis of similarity (ANOSIM) to test the significant differences between altitude gradients. NMDS and ANOSIM tests were calculated using the "metaMDS" and "anosim" functions. The "indicspecies" package was used to investigate the association between species and altitude by identifying the indicator soil surface-dwelling Oribatida species of each altitude (Cáceres & Legendre, 2009;Dufrêne & Legendre, 1997).
After excluding the rare species (abundance <1%), variation partitioning analysis (VPA) was performed to explore the effect of environmental factors (including TC, TN, TP, TH, TS, and QUA), climate factors (including MAT and MAP), and an elevation factor (ALT) on variations of soil surface-dwelling Oribatida community. The VPA was calculated with the "varpart" function in the "vegan" package, then used the "rda" function to test the significance of VPA results.
We tested the correlation of spatial factors (including geographic distance and elevation distance) and soil surface-dwelling Oribatida community using the Mantel analysis (Legendre & Legendre, 2012).
The Mantel test was calculated with the "mantel" function in the "vegan" package. The geographic distance between each altitude was calculated with the "distm" function in the "geosphere" package (Hijmans et al., 2019). In addition, we used the ArcGIS 10.2 software to plot the study area and sampling sites. All data analyses were performed in R 4.2.0 (R Core Team, 2022).

| Alpha diversity
A total of 34 species of soil surface-dwelling Oribatida were identified along the altitudinal gradient (  Figure 2a,b). Shannon-Wiener and Simpson indexes were significantly higher at 800 m than at other altitudes (Kruskal-Wallis test, p < .001; Figure 2c,d). The density and alpha diversity of soil surface-dwelling Oribatida community showed a significantly negative correlation with altitude ( Figure 3).

| Beta diversity
Beta diversity of soil surface-dwelling Oribatida was significantly lower at 2300 m than at other altitudes (p < .001, Figure 4). The partition in beta diversity revealed that differences in species composition were generally induced by species turnover and not nestedness.
The 600 m, 800 m, 1600 m, and 2000 m altitudes were significantly influenced by the species turnover component, and the 2300 m altitude was significantly influenced by the species nestedness component.

| Community assemblages
The non-metric multidimensional scaling (NMDS) ordination iden- showed that Oribatida assemblages significantly differed along the altitudinal gradient (Table 3). Indicator analysis revealed representative taxa for each altitude gradient (Table 4)

| Variation of community
Nine factors were selected as the drivers of variation in the Oribatida community structure (Table 1). Venn diagram showed that elevation factor (altitude) and climate factors (including MAT and MAP) had an extremely significant correlation on the soil surface-dwelling Oribatida community (Figure 6, p < .01). The most substantial contribution was attributed to the elevation factor (22%).
Soil surface-dwelling Oribatida community dissimilarity had a positive correlation with geographic and elevation distance, explaining 27.08% and 25.63% of the variations, respectively (Figure 7). The variation of soil surface-dwelling Oribatida community structure between 600 and 800 m and between 600 and 2300 m had a significant correlation with geographic distance (Table 5). The variation in soil surface-dwelling Oribatida community structure between 1600, 2000, and 2300 m significantly correlated with elevation distance.

| DISCUSS ION
Our study assessed the patterns and drivers of soil surface-dwelling Oribatida diversity along an altitudinal gradient on the Changbai Mountain. We found that the distribution pattern of alpha diversity and beta diversity of soil surface-dwelling Oribatida decreased with the rising altitude. Species turnover was the essential component in driving beta diversity. NMDS analysis showed that the community structure of soil surface-dwelling Oribatida assemblages varied significantly along the altitudinal gradient. The community variation of soil surface-dwelling Oribatida was affected by altitude and climate (temperature and precipitation). Geographical and elevation distances were the primary drivers of variation in the distribution patterns of soil surface-dwelling Oribatida communities.

| Alpha diversity
The alpha diversity of soil surface-dwelling Oribatida assemblages showed a decreasing trend with the rising altitude. However, another peak of soil surface-dwelling Oribatida density was recorded at 2000 m rather than middle altitude. In our study, the sampling sites in birch forests are located at a middle-high altitude of approximately 2000 m. A balance has been established between increased precipitation and decreased temperature, and relatively high-quality litter per unit area had provided a diverse microhabitat for soil surface-dwelling Oribatida (Hasegawa et al., 2006;Xu et al., 2020).

Such a suitable habitat makes three species of soil surface-dwelling
Oribatida (E. borealis, P. peltifer, T. berlesei) reach an extremely high density, resulting in another density peak at 2000 m. This distribution pattern can also relate to edge effects, as the forest line is the most prominent biogeographic boundary of the montane ecosystem (Becker et al., 2007). The sampling area of 2000 m is a patchy transition zone consisting of several transitional vegetation zones and their components. It is near the timberline between the Changbai Mountain forest zone and the alpine tundra. Such habitat boundary often shows a high plant diversity, biomass, species dispersal, and soil nutrient resources (Reinmann & Hutyra, 2017;Ruwanza, 2018;Shen et al., 2021). For instance, compared to other vegetation zones, the microbial biomass in the Changbai Mountain forest transition zone was higher (Liu et al., 2019).
The environmental filters generally become more robust with an increased altitude (Hoiss et al., 2012;Webb et al., 2002;. The soil surface-dwelling Oribatida were dominated by a combination of significantly highest nestedness at 2300 m due to lack of resources, harsh climatic conditions, and topographic barriers (Castro et al., 2019;Melo et al., 2009).
In our study, the 2300 m (tundra zone) altitude was dominated by two soil surface-dwelling Oribatida species, P. peltifer and T. berlesei, which accounted for 89.72% of the total community density, associated with the loss of species occurred in the high altitude (Willig & Presley, 2016). P. peltifer and T. berlesei were not the dominant species at 600, 800, and 1600 m because a common species widely distributed at high altitudes could usually only be active in some small extreme patches in a lower altitude (Otto & Svensson, 1982).

F I G U R E 3
Multiple linear regressions between altitude and density, richness, Shannon-Wiener index, and Simpson index of soil surfacedwelling Oribatida along an altitudinal gradient on the Changbai Mountain, Northeast China.

F I G U R E 4
Beta diversity (a) and its turnover (b) and nestedness (c) components based on Sørensen dissimilarity between five altitudes. Different letters above boxes indicate significant differences between altitudes (p < .05, Kruskal-Wallis test).
In our study, we observed an average density of 3286 ind./m 2 for soil surface-dwelling Oribatida, which is considerably lower than the density of soil-dwelling Oribatida reported by Hasegawa et al. (2006) and Illig et al. (2010), ranging from 15,000 to 25,000 ind./m 2 . The sampling tool had a large impact on experimental results (Jing et al., 2005). We compared with Liu et al. (2023) and discovered that a proportion of species of Oribatida were exclusively collected using the Vortis arthropod suction sampler, such as Scapheremaeus polysetosus, Archoplophora rostralis, Micreremus brevipes and Oribotritia gigas, and several Oribatida were exclusively collected using Tullgren extractors, such as Parhypochthonius aphidinus, Liochthonius sellnicki and Synchthonius elegans. Therefore, combining with different sampling tools is a better option in fieldwork, for a more accurate and comprehensive collection of sample data.

| Beta diversity
From 800 to 2300 m, the overall beta diversity of soil surfacedwelling Oribatida decreased gradually. A significant diversity variation of soil surface-dwelling Oribatida occurred at 2300 m, which also supported our first Hypothesis. The overall beta diversity was significantly higher at 600, 800, and 1600 m than at other altitudes. In such cases, even a slight increase in the area could incorporate many new species of soil surface-dwelling Oribatida (Jankowski et al., 2009). The species turnover component was the primary process of soil surface-dwelling Oribatida community structuring, which supported our second Hypothesis. This pattern of variation dominated by species turnover reflects a deterministic process of community construction, consistent with the study by Bishop et al. (2015).
Some species have a remarkable dispersal capacity and can tolerate the frequent extreme climates and low availability of resources in high altitudes (Hasegawa et al., 2006;Illig et al., 2010;Otto & Svensson, 1982;Valencia-Cuevas et al., 2020).

F I G U R E 6
The contribution of abiotic parameters to soil surface-dwelling Oribatida abundance along an altitudinal gradient on the Changbai Mountain, Northeast China (ALT, altitude; ENV, environmental factors including total carbon, total nitrogen, total hydrogen, total sulfur, total phosphorus, and the litter quality of per unit area; MAP, mean annual precipitation; MAT, mean annual temperature), **p < .01.
acuta, Lauroppia neerlandica, among others) could establish at all altitudes. In contrast, other species with narrower ecological ranges can only present at lower altitudes. This pattern reflected the substantial turnover of soil surface-dwelling Oribatida (Marian et al., 2018).
Because of the adaptive divergence of species along the altitudinal gradient, species living in low altitudes lack the adaption mechanism for middle-high altitudes (Keller et al., 2013;Samson et al., 1997).

| Variation of community
In our study, spatial filtering strongly correlated with the soil surface-dwelling Oribatida community, consistent with our third Hypothesis. The elevation is a vital compound variable in montane ecosystems (Mumladze et al., 2017). We found that the elevation had an extremely significant effect on the community structure of soil surface-dwelling Oribatida. However, other studies found that the elevation had no biological significance, and the variation of the microarthropod community along the altitudinal gradient was associated with either biotic or abiotic factors, such as plant community, rather than a response to altitude (Jiang et al., 2015;Nieto Peñalver et al., 2017;Valencia-Cuevas et al., 2020).
Microarthropod communities are significantly associated with climatic factors, composition and quality of litter, soil characteristics including pH, and resource availability (Shen et al., 2013;Winkler et al., 2018;. Our study revealed that soil surface-dwelling Oribatida were not significantly influenced by environmental factors, except for temperature and precipitation. Temperature is an indirect measure of primary productivity and a determinant of arthropod community composition, representing the altitudinal limit of species distribution in the montane eco- other ecosystems have also demonstrated the effect of precipitation on Oribatida (Hense, 2016;Lehmitz et al., 2020). Our study can support the third Hypothesis that temperature and precipitation could be the primary drivers of soil surface-dwelling Oribatida variation along the altitudinal gradient. Microarthropod diversity is closely associated with carbon content and carbon release (Fujii et al., 2018;Liu et al., 2019), and the nitrogen content can impact the feeding behavior of microarthropods by limiting microbial communities (Devetter et al., 2017). However, the total carbon and total nitrogen content did not significantly affect the abundance of microarthropods in our study.

| CON CLUS IONS
In conclusion, our results revealed a dynamic pattern of soil surfacedwelling Oribatida diversity and community structure in montane ecosystems along altitudinal gradients. The abundance of soil surface-dwelling Oribatida can be influenced by elevation and climate factors. Spatial factors, including geographic and elevation distance, were the primary factors associated with the variations of soil surface-dwelling Oribatida community in montane ecosystems.
Our study provides evidence for the biodiversity analyzing of soil surface-dwelling Oribatida in montane ecosystems along an altitudinal gradient.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declared that they have no conflicts of interest in this work.